
function main_jk()


% sets some parameters
k=2.25;  % sets bandwidth proportionality
g1=[(0.399999:.01:0.499999)';(0.500001:.01:0.650001)'];
g2=[(0.399999:.01:0.499999)';(0.500001:.01:0.650001)'];

% loads data explorted from stata
by=importdata('D:\research\connan\compplan\data\matlab\dd_lccpresence.txt');
bx=importdata('D:\research\connan\compplan\data\matlab\dd_lcc_top2.txt');
bn=size(by,1);

o_id=importdata('D:\research\connan\compplan\data\matlab\dd_lcc_origin_id.txt');
d_id=importdata('D:\research\connan\compplan\data\matlab\dd_lcc_dest_id.txt');


% calculates jack-knife standard errors
ps_yhat=zeros(27,27,bn);
ps_g=zeros(27,2);

for i=1:bn
    
    % leaves out one block at a time
    ind=(o_id~=o_id(i) & d_id~=d_id(i));
    psy=by(ind);
    psx=bx(ind,:);
    clear ind;
    
    % chops up data exported from stata
    [y1,x1,n1,y2,x2,n2,y3,x3,n3,y4,x4,n4] = gen_data(psy,psx);
    
    y_hat = local_linear_rdd(y1,x1,n1,y2,x2,n2,y3,x3,n3,y4,x4,n4,g1,g2,k);
    clear y1 x1 n1 y2 x2 n2 y3 x3 n3 y4 x4 n4;
    
    ps_yhat(:,:,i)=y_hat;
    ps_g(:,1)=g1;
    ps_g(:,2)=g2;
    clear y_hat;
    
end


std_error=std(ps_yhat,[],3);
std_error1=std(diff(ps_yhat,1,1),[],3);
std_error2=std(diff(ps_yhat,1,2),[],3);

save stderr std_error std_error1 std_error2 ps_yhat ps_g;


end
